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ABSTRACT 

We adapt the L-Galaxies semi-analytic model to follow the star formation histories 
(SFH) of galaxies - by which we mean a record of the formation time and metallicities of the 
stars that are present in each galaxy at a given time. We use these to construct stellar spectra in 
post-processing, which offers large efficiency savings and allows user-defined spectral bands 
and dust models to be applied to data stored in the Millennium data repository. 

We contrast model SFHs from the Millennium Simulation with observed ones from the 
VESPA algorithm as applied to the SDSS-7 catalogue. The overall agreement is good, with 
both simulated and SDSS galaxies showing a steeper SFH with increased stellar mass. The 
SFHs of blue and red galaxies, however, show poor agreement between data and simulations, 
which may indicate that the termination of star formation is too abrupt in the models. 

The mean star-formation rate (SFR) of model galaxies is well-defined and is accurately 
modelled by a double power law at all redshifts: SFR ex 1/(a; -1 ' 39 + at 1 ’ 33 ), where x = (t a — 
t)/'■>.{) Gyr, t, is the age of the stars and t H is the lookback time to the onset of galaxy formation; 
above a redshift of unity, this is well approximated by a gamma function: SFRoc x 15 e~ x , 
where x = (t a — t) /2.0 Gyr. Individual galaxies, however, show a wide dispersion about this 
mean. When split by mass, the SFR peaks earlier for high-mass galaxies than for lower-mass 
ones, and we interpret this downsizing as a mass-dependence in the evolution of the quenched 
fraction: the SFHs of star-forming galaxies show only a weak mass dependence. 
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1 INTRODUCTION 

Understanding the astrophysics behind the formation and evolution 
of galaxies is an important goal in modern astronomy. One of the 
most fundamental probes of that physics is the star formation rate 
as a function of cosmic time. In this paper we contrast predicted 
and observed star formation histories (hereafter, SFHs) of galaxies, 
and explore the expected range of SFHs at high redshift. 

Two main observational approaches are used to infer the SFH 
of galaxies. One can look at the instantaneous star formation rate 
as a function of cosmic time (for an overview see, e.g. Kennicutt 
1998; Calzetti 1999), or one can deduce the history from the fossil 
record of current-day galaxies. The two techniques in fact measure 
slightly different things, with the relationship between the two de¬ 
pending upon the merging history of the galaxies. 

* E-mail: P.A.Thomas@sussex.ac.uk 


This paper focuses upon the archaeological approach. We have 
extended the L-GALAXIES semi-analytic model (hereafter simply 
L-GALAXIES) to keep a record of the SFHs of individual galaxies, 
with a bin-size that increases with increasing age for the stars. The 
resulting data have been made available as part of the public data 
release (DR) of the Millennium Simulation 1 (Lemson & The Virgo 
Consortium 2006) that accompanies the latest implementation of 
L-GALAXIES (Henriques et al. 2014). 

We note that the SFH as defined in this paper (the history of 
star formation of all the stars that end up in the galaxy at some par¬ 
ticular time) is related to, but distinct from, the variation in mass of 
a galaxy over time. The latter follows only the history of the main 
component of the galaxy (along the “main branch” of the merger 
tree) and has been investigated for the Millennium Simulation by 


1 http://gavo.mpa-garching.mpg.de/MyMillennium/ 
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Cohn & van de Voort (2015). The difference between the two re¬ 
flects the merger history of galaxies. 

The term SFH is often loosely used in papers without being 
defined. Observationally, the only direct measure of SFHs corre¬ 
sponds to that described in this paper, i.e. the distribution of for¬ 
mation times of all the stars that make up the galaxy. 2 Measures of 
the SFFI of the main galactic component can only be inferred sta¬ 
tistically by observing populations of galaxies at different redshifts 
and making some assumptions about merger rates as a function of 
stellar mass and environment. In principle, the former method is 
much cleaner, as it is free from these model assumptions; however, 
in practice, the inversion of the stellar spectra is highly degenerate 
and sensitive to the input stellar population synthesis models, and 
can lead to implausible results if some model constraints are not 
imposed. 

As well as enabling comparison to observations, the introduc¬ 
tion of SFFls into L-GALAXIES allows reconstruction of galaxy 
magnitudes using arbitrary stellar population synthesis and dust 
models, in any band, in post-processing, and we investigate the ac¬ 
curacy of this approach. In addition, having SFHs is a prerequisite 
for a correct, time-resolved treatment of galactic chemical enrich¬ 
ment, as described in Yates et al. (2013). 

The star formation history of a galaxy, its chemical evolution, 
and its current dust content, can in principle be fully recovered with 
high enough quality observational data, suitable modelling of the 
spectral energy distribution of stellar populations and dust extinc¬ 
tion, and appropriate parametrization. Several algorithms have been 
developed over the last decade that attempt to do the above in the 
most robust way (see e.g. MOPED by Heavens et al. 2004; Pan- 
ter et al. 2007, STECMAP by Ocvirk et al. 2006, STARLIGHT by 
Cid Fernandes et al. 2004, 2005 or ULYSS by Koleva et al. 2009). 
In this paper we focus on the results obtained by VESPA (Tojeiro 
et al. 2007), a full spectral fitting code that was applied to over 
800,000 Sloan Digital Sky Survey DR7 galaxies (SDSS Collabora¬ 
tion 2009). The resulting data base of individual SFHs, metallicity 
histories and dust content is publicly available and described in To¬ 
jeiro et al. (2009, hereafter TWH09). The wide range of galaxies 
in SDSS DR7 and the time resolution of the SFHs published in the 
data base make it ideal for a detailed comparison with our model 
predictions. 

At higher redshift, observational data are scarce. Nevertheless, 
we make predictions from our models of how the star formation rate 
of galaxies evolves with redshift, and we show how downsizing 
arises from a mass-dependence of the rate at which galaxies are 
quenched. 

In outline, the main aims of this paper are: 

• to briefly overview the L-GALAXIES and VESPA algorithms 
- Sections 2.1 and 2.2; 

• to describe the L-GALAXIES SFH binning method - Sec¬ 
tion 2.3; 

• to test how well post-processing can reconstruct magnitudes - 
Section 3; 

• to compare model galaxies from L-GALAXIES with results 
from the VESPA catalogue - Section 4; 

• to investigate the variety of SFHs that we find in our model 
galaxies - Section 5; 

• to provide a summary of our key results - Section 6. 

2 We note that this would better be described by the phrase “stellar age 
distribution” rather than star formation history, but the latter phrase pre¬ 
dominates in the literature. 


2 METHODS 

2.1 L-Galaxies 

In this work we use the latest version of the Munich semi-analytic 
(SA) model, L-GALAXIES, as described in Henriques et al. (2014). 
This gives a good fit to, amongst other things, the observed evolu¬ 
tion of the mass and luminosity functions of galaxies, the fraction 
of quenched galaxies, the star formation versus stellar mass relation 
(at least at z < 2), the Tully-Fisher relation, metallicities, black 
hole masses, etc. We refer to this standard model as HWT14. The 
improvements obtained in terms of the evolution of the abundance 
and red fractions of galaxies as a function of stellar mass in this 
model result from: a supernova feedback model in which ejected 
gas is allowed to fall back on to the galaxy on a time-scale that 
scales inversely with halo virial mass (Henriques et al. 2013); and 
a lower star formation threshold and weaker environmental effects 
both reducing the suppression of star formation in dwarf galaxies. 

L-GALAXIES, as with other SA models, follows the growth 
of galaxies within the framework of a merger tree of dark mat¬ 
ter haloes. We construct this tree from the Millennium Simulation 
(Springel et al. 2005), scaled using the method of Angulo & White 
(2010) to the Planck cosmology (Planck Collaboration XVI 2014): 
H m = 0.315, Ha = 0.683, H b = 0.0488, n s = 0.958, erg = 
0.826, h = 0.673. This then gives a box size of 480.3 /i -1 Mpc 
and a particle mass of 9.61 x 10 8 /i _1 Mq. The merger tree is con¬ 
structed from 58 snapshots, 3 each of which is subdivided into 20 
integration timesteps. The snapshots are unevenly spaced, such that 
the time resolution is higher at high redshift, but a typical timestep 
is of the order of 1-2 x 10' yr. We note that the SA model has been 
implemented on both the higher resolution Millennium-II (Boylan- 
Kolchin et al. 2009) and larger volume Millennium-XXL simula¬ 
tions (Angulo et al. 2014) although we do not make use of those 
simulations in the current work. 

Prior to HWT14, galaxy luminosities were determined by ac¬ 
cumulating flux in different spectral bands throughout the time- 
evolution of each galaxy. In the new model those fluxes are recov¬ 
ered to high accuracy in post-processing, by simply recording the 
star formation and metallicity history in a relatively small number 
of time-bins. The new method is introduced following the work out¬ 
lined and presented in the current paper. The DR that accompanies 
this series of papers records the SFHs, allowing the user flexibility 
to define their own bands and dust models. 

2.2 VESPA 

The spectrum of galaxy, in the absence of dust, can be described 
as the linear superposition of the spectra of the stellar populations 
of different ages and metallicities that exist in the galaxy. The de- 
convolution of a galaxy’s spectrum into a star formation and metal¬ 
licity history is in principle trivial, but complicated by noisy or in¬ 
complete data and limitations in the modelling. Ocvirk et al. (2006) 
showed how the problem quickly becomes ill conditioned as noise 
increases in the data, and that the risk of overparametrizing a galaxy 
is high. VESPA takes into account the noise and data quality of each 
individual galaxy and uses an algebraic approach to estimate how 
many linearly independent components one can extract from each 
observed spectrum, thereby avoiding fitting the noise rather than 
the signal - see Tojeiro et al. (2007) for details. 

3 Five of the original Millennium Simulation snapshots lie at z < 0 after 
scaling. 
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Figure 1 . The evolution of the SFH bins for N m ax = 2. The at-axis rep¬ 
resents successive timesteps, starting from high redshift and moving to¬ 
wards the present day. The y-axis represents bins of stellar age, starting 
with stars that are newly-created at that time and looking back into the past. 
The numbers within each bin represent the number of timesteps that have 
been merged together to produce that bin. The columns that are bracketed 
together show different arrangements of the data at a single cosmic time. 
The shaded, red stacks represent transient structures in which some bins 
merge together to produce the black stacks to their right. 

VESPA recovers the SFH of a galaxy in 3 to 16 age bins 
(depending on the quality of the spectra), logarithmically spaced 
between 0.002 Gyr and the age of the Universe. For each age 
bin, VESPA returns the total mass formed within the bin and a 
mass-weighted metallicity of the bin, together with an estimate of 
the dust content of the galaxy. As we always compare our model 
predictions to the mean SFH of large ensembles of galaxies, we 
choose to use the fully-resolved SFHs published in the data base of 
TWH09. Whereas we expect these to be dominated by the noise on 
each individual galaxy, the mean over a large ensemble has been 
shown to be robust (Panter et al. 2003). In this paper, we will com¬ 
pare to mean SFHs for different galaxy populations, as described 
in the text. 

2.3 The SFH binning algorithm 

As mentioned in Section 2.1, the Millennium merger trees are con¬ 
structed from 58 snapshots, each of which is separated into 20 inte¬ 
gration timesteps. To follow the history of star formation, we intro¬ 
duce extra arrays to carry information on the mass and metallicity 
of stars in each component of the galaxy (disc, bulge, intracluster 
mass) as a function of cosmic time. To save the information over 
all 1160 timesteps would consume too much memory and is unnec¬ 
essary. Instead, we wish to use a high resolution for the recent past 
(when the stellar population is rapidly evolving) and a lower one 
at more distant times. To do so, we adopt the following procedure 
(see Fig. 1). 

Starting at high redshift, on each timestep a new bin is cre¬ 
ated to hold the SFH information. Whenever the number of bins 
of a particular resolution exceeds N max (where in the diagram, 
for the purposes of illustration, N max = 2), then the two oldest 
bins are merged together to form a new bin of twice the size - this 
may result in a cascade of mergers at successively higher levels 
(in the figure, these mergers are represented by red columns joined 
by braces to the merger product). In this way, the number of bins at 


Table 1 . For the Millennium Simulation, using 1160 timesteps, this table 
shows, for different choices of V max : IVtot - the maximum number of SFH 
bins required; N z= o — the number of populated bins at z = 0: Af m j n /yr - 
the minimum bin-size in years at z = 0; Af max /yr - the maximum bin-size 
in years at z = 0. 


N max 

IVtot 

O 

II 

£ 

Af min /yr 

A^max/yr 

i 

11 

1 

6.0 x 10 7 

1.1 x 10 10 

2 

19 

16 

1.5 x 10 7 

2.1 x 10 9 

3 

27 

23 

1.5 x 10 7 

1.6 x 10 9 

4 

34 

31 

1.5 x 10 7 

5.6 x 10 8 


each size grows from 1 to (V max , then oscillates between JV max and 
IV max — 1. The total number of bins required does not exceed the 
smallest integer greater than AUax log 2 (W s tep/A r max + 1), where 
IVstep is the number of timesteps. Table 1 shows, for the Millennium 
Simulation, using 20 steps within each of 58 snapshots, the maxi¬ 
mum number of bins required, and, at z = 0, the actual number 
of SFH bins and their minimum and maximum size in years. 4 Note 
that all choices of IV max > 2 have the same minimum bin-size, 
equal to that of the original timesteps; what differs is the number of 
bins that are resolved at that highest resolution. 

We have investigated the sensitivity of our results to the num¬ 
ber of bins and conclude that iV max = 2 gives the best balance 
between data-size and accuracy: that is the value used in the Mil¬ 
lennium data base and, unless mentioned otherwise, in the results 
presented below. 


3 POST-PROCESSING OF MAGNITUDES 

In most SA models, and in F-GALAXIES prior to this work, galaxy 
luminosities are computed by adding the flux in different bands 
throughout the time-evolution of each galaxy. This calculation gen¬ 
erally requires interpolating between values in large stellar popula¬ 
tion synthesis tables and represents a large fraction, and in some 
cases the majority, of the computational time for the entire galaxy 
formation model. The problem is aggravated as different types of 
magnitudes (dust corrected, observer-frame) for additional compo¬ 
nents (e.g. disc, bulge, intra-cluster light) are included. 

These difficulties can, in principle, be circumvented by stor¬ 
ing the star formation and metallicity histories for different compo¬ 
nents of the galaxies, and using them to compute emission in post¬ 
processing. Ideally, this history would be stored for all the interme¬ 
diate steps between output snapshots for which galaxy properties 
are computed. However, memory constraints make this infeasible. 
For our current set up, for example, it would require storage of up 
to 2320 values for each galaxy component (58 snapshots, 20 inter¬ 
mediate steps per snapshot, star formation and metallicity). 

An alternative, tested in this paper, is to store the histories in 
bins that grow in size for older populations, as described in Sec¬ 
tion 2.3. Since the emission properties of populations vary on sig¬ 
nificantly longer time scales for old populations this can in princi¬ 
ple allow us to maintain accuracy. To validate the method we com¬ 
pare the theoretical emission from galaxies computed both with full 
resolution on-the-fly and from star formation and metallicity his¬ 
tory bins that are merged together for older populations. We assume 

4 These numbers are unchanged for the 63 snapshots required to use the 
Millennium Simulation with the original WMAP -1 cosmology. 
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Table 2. This table shows the root-mean-square difference between magnitudes calculated on-the-fly during the running of the code and those calculated in 
post-processing. 


2 = 0 


Without dust 



With dust 


N max — 1 

A^max — 2 

A'max — 4 

Nmax — 1 

Afmax — 2 

N max — 4 

GALEX FUV 

3.02 

0.29 

0.16 

2.53 

0.38 

0.29 

GALEX NUV 

1.79 

0.05 

0.03 

1.47 

0.16 

0.16 

SDSS u 

0.76 

0.03 

0.01 

0.66 

0.06 

0.05 

SDSS g 

0.50 

0.02 

0.00 

0.46 

0.03 

0.02 

SDSS 2 

0.36 

0.02 

0.00 

0.33 

0.02 

0.02 

VISTA J 

0.35 

0.02 

0.00 

0.32 

0.02 

0.02 

VISTA K s 

0.33 

0.01 

0.01 

0.31 

0.02 

0.02 

IRAC 3.6 

0.38 

0.01 

0.00 

0.36 

0.03 

0.02 



without dust 



with dust 


2-2 








Afmax — 1 

N max — 2 

AW = 4 

A^max = 1 

N max — 2 

AT max = 4 

GALEX FUV 

1.74 

0.04 

0.01 

1.14 

0.27 

0.26 

GALEX NUV 

1.39 

0.02 

0.00 

0.95 

0.21 

0.20 

SDSS u 

0.88 

0.02 

0.00 

0.63 

0.13 

0.12 

SDSS g 

0.64 

0.01 

0.00 

0.47 

0.09 

0.09 

SDSS 2 

0.39 

0.01 

0.00 

0.33 

0.04 

0.04 

VISTA J 

0.33 

0.01 

0.00 

0.30 

0.02 

0.02 

VISTA K s 

0.32 

0.01 

0.00 

0.30 

0.01 

0.01 

IRAC 3.6/im 

0.31 

0.01 

0.00 

0.28 

0.02 

0.02 




- 2 - 1012 - 2-1012 - 2 - 1012 - 2-1012 
^ ® 9 tobies — M C3 9 post_process ^ ^ 9 tobies — M O g p0 st_process ^ ® 9 tobies — ^ ® 9 post-process 9 tobies - ^^9 post-process 

Figure 2. The difference between photometric properties calculated using full resolution and binned SFHs at z = 0 (left) and z = 2 (right). The left-hand 
column in each set of panels correspond to intrinsic magnitudes while the right takes into account dust extinction. From top to bottom the resolution of the 
binning is increased: N max = 1 (top), 2 (middle) and 4 (bottom). 
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that star formation occurs at a time corresponding to the mid-point 
of each SFH bin. To spread star formation out over the time-bin 
would be equivalent to using a larger number of timesteps (which 
we have also tested) and makes little difference except in the UV. 

Fig. 2 shows the difference between photometric properties 
calculated using full resolution and binned SFHs for 2 = 0 and 
2 = 2. In both sets of panels the left-hand column corresponds 
to intrinsic magnitudes while the right-hand column takes into ac¬ 
count dust extinction. From top to bottom the resolution of the bin¬ 
ning is increased: IV max = 1, 2 and 4. While large differences 
between the two methods are seen for the lowest resolution, the 
figures show that good convergence is achieved for IV max = 2 or 
more. 

Table 2 shows the root-mean-square (rms) difference between 
magnitudes calculated on-the-fly during the running of the code 
(i.e. using the finest possible time resolution) and those calculated 
in post-processing, using different numbers of SFH bins. Quantita¬ 
tively, for (V max = 2, at 2 = 0, the rms difference between the two 
methods is less than 0.05 for all intrinsic magnitudes except the 
far-UV, for which it is approximately 0.29. 5 At 2 = 2, the mean 
difference is less than 0.04 in all bands. 

The increased accuracy in the far-UV at high redshift results 
from the higher accuracy in post-processing at a time in which 
fewer bins were merged. Emission in this part of the spectrum 
is dominated by extremely young populations for which even a 
slightly different formation time results in a large variation in pre¬ 
dicted flux. The two methods thus differ in detail, but have the same 
statistical properties. If one is interested in the precise UV flux of 
a particular galaxy then that can be recovered by using finer time- 
bins. We have checked that, keeping 7V max = 2 but using a finer 
timestep (which adds very few SFH bins) improves the agreement 
between the two methods of calculating fluxes. 6 

At low redshift, it can be seen that there is a residual error 
in the calculated intrinsic far-UV flux even for iV max = 4. This 
arises from stars of age about 1 Gyr (i.e. the TP-AGB population), 
for which the tabulated UV fluxes in the Maraston (2005) popu¬ 
lation synthesis tables show a large jump in luminosity between 
the two lowest-metallicity bins. Merging galaxies that contain stars 
of differing metallicity can therefore lead to large changes in flux. 
Such merging can occur for any choice of A r max and is hence a 
fundamental (albeit very minor) limitation of the SFH magnitude- 
reconstruction method. 7 

The addition of dust significantly degrades the agreement be¬ 
tween the two reconstruction methods in all bands, with the far- 
and near-UV being most affected. At 2 = 2 the rms differences 
between dust-corrected magnitudes are approximately 0.20 for far- 
and near-UV, 0.10 for u and g, 0.04 for 2 and J, and 0.02 for K s 
and irac3.6/im- 

In the current version of L-GALAXIES a two-component dust 
model applies extinction separately from the diffuse interstellar 
medium and from molecular birth clouds (see Section 1.14 in the 
supplementary material of Henriques et al. 2014 for details). The 
large differences seen for dust corrected magnitudes are mostly 

5 These values are reduced to 0.03 and 0.16, respectively, when V max = 4 
is used. 

6 We do not use a finer time resolution as a default as this agreement is 
illusory — the underlying merger tree is not capturing the dynamics on that 
short a timescale. 

' The limitation could be overcome by keeping SFHs for several different 
metallicity bins, but this moves away from the spirit of the method and we 
do not think that the gain justifies the extra storage cost. 


caused by the latter. The calculation of this optical depth includes 
a random gaussian term that leads to differences in the amount ex¬ 
tinction assumed for each individual galaxy when a different num¬ 
ber of timesteps are used. 

The method successfully tested in this section is adopted in the 
recent major release of the Munich model, HWT14. By computing 
emission properties in post-processing, the memory consumption 
of the code is no longer dependent on the number of photomet¬ 
ric bands. Moreover, the method allows emission properties to be 
computed after the model is completed using any stellar popula¬ 
tions synthesis code and for the filters used by any observational 
instrument. To show the potential of the new method, the new ma¬ 
jor release already includes emission in 20 bands in the snapshot 
catalogues, and in 40 bands and for two different stellar popula¬ 
tions for the lightcones, all calculated in post-processing. 


4 COMPARISON OF VESPA AND L-Galaxies 

For reasons described in Appendix A, we use the version of the 
VESPA catalogue that was created using the population synthesis 
model of Maraston (2005) with a one-component dust model. 

Note that the VESPA data and the SA models produce output 
with very different binning. The time resolution in the observations 
is necessarily very coarse at high redshift, whereas there is no such 
restriction in the models. In Section 4.1 below, we re-bin the model 
predictions to match those of VESPA; throughout the rest of the 
paper, we will keep the actual binning returned by the models so as 
to allow a clearer understanding of the growth of galaxies at high 
redshift. 

4.1 The main galaxy sample 

The main SDSS galaxy sample covers a redshift range of 0 < 2 < 
0.35. In Fig. 3 we show the mean SFH for all galaxies using the 
maximum VESPA resolution of 16 bins. For most galaxies, the data 
quality is not good enough to independently measure masses in all 
16 bins and so the VESPA algorithm will return solutions on bins 
of varying and lower-resolution width, as described in TWH09. 
The assumed star formation rate (hereafter, SFR) within each bin 
depends on its width (it is constant in high-resolution bins, and 
exponentially decaying in low-resolution bins). Choosing a SFR 
within a bin is an unavoidable part of the process of parametrizing 
a galaxy. We have however checked that our conclusions remain 
unchanged if we: (i) use only galaxies with high-resolution bins, 
and (ii) use a constant SFR in wide bins. 

For each VESPA galaxy, we select the model galaxy that most 
closely matches it in mass and redshift, then use this to construct a 
mean SFH. The result is shown in Fig. 3 along with the predictions 
from two versions of L-GALAXIES: GWBll-Guo et al. (2011, ma¬ 
genta diamonds and dot-dashed line) and HWT14-Henriques et al. 
(2014, green squares and dashed lines). 

In this and subsequent plots, A/ is defined as the fraction of 
stars (i.e. the specific stellar mass) within each bin. Unless men¬ 
tioned otherwise, A/ is calculated separately for each galaxy and 
then averaged (i.e. weighted by galaxy number rather than stellar 
mass). 

The VESPA results are not as smooth as the SA model data. 
We expect the averaged rest-frame SFH of a varied ensemble of 
galaxies to be devoid of significant structure, and to appear smooth 
as seen on the model galaxies. At late times, especially, the SFR 
should be approximately constant, and so one would expect the 
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Figure 3. The average SFH from VESPA (black lines), the GWB11 model 
(magenta, dash-dotted lines) and the HWT14 model (dashed, green lines) 
within the indicated mass range. To guide the eye, the blue, dotted line is the 
same in each panel and has slope unity, corresponding to a constant SFR. 

SFH bins to run parallel to the dotted line, of slope equal to unity. 
The features seen on the data are heavily dependent on the mod¬ 
elling (see e.g. TWH09, Tojeiro et al. 2013, Appendix A), and 
therefore are the likely result of limitations of the stellar popula¬ 
tion synthesis and dust modelling. 

The SA models show a turn-down in SFR in the oldest stel¬ 
lar age bin, corresponding to the onset of star formation. No such 
feature is seen in the VESPA results, most probably because the 
spectral signatures are too weak to be detected and so the recon¬ 
struction method imposes a declining SFR even at these earliest 
times. 

Other than that, at ages above 1 Gyr the GWB11 model 
seems to provide a reasonable match the observations, whereas for 
younger stars, the HWT14 model is a better fit. To draw more 
definitive conclusions about which is the preferred model, one 
would have to look in much more detail at the reconstruction bi¬ 
ases that may be present in the VESPA method when applied to 
imperfect data, and that will be the subject of future work. 

From here on we will rebin the VESPA solutions to five bins 
in age, to reduce the scatter and to average over these features that 
we know to be unphysical. 

4.2 Mass selection 

Fig. 4 shows the SFHs broken down by stellar mass. The VESPA 
reconstruction gives a slope that is too steep for high-mass galax¬ 
ies and too shallow for low-mass galaxies, when compared to the 
expected constant SFR at recent times. Nevertheless, it can be seen 
that both the VESPA galaxies, and those in the SA model, form 
stars earlier in higher-mass galaxies and have a correspondingly 
lower SFR at late times. 

The variation of the model SFHs with mass is explored further 
in Section 5.3, below. 

4.3 Colour selection 

Next, we look at the distinction between red and blue galaxies by 
selecting according to u — r colour. As HW14 showed, in spite of 
reproducing the observed galaxy colour bimodality, L-GALAXIES 



0.01 0.1 1 10 
Stellar age / Gyr 


Figure 4. The average SFH from VESPA (black lines), the GWB 11 model 
(magenta, dash-dotted lines) and the HWT14 model (dashed, green lines) 
within two different stellar mass bins, as shown. To guide the eye, the blue, 
dotted line is the same in each panel and has slope unity, corresponding to 
a constant SFR. 

does not reproduce the exact colour distributions seen in SDSS. 
Therefore, applying the same colour cuts in the data and simula¬ 
tion would result in picking out intrinsically different galaxy pop¬ 
ulations. Instead, we select the 10% bluest and reddest galaxies, 
according to u — r colour, in both the VESPA and L-GALAXIES 
samples. The resulting SFHs are show in Fig. 5 for galaxies of mass 
10 10 <M/Mq< 10 10 ' 5 at low redshift, z < 0.07. 

VESPA produces similar SFHs for both blue and red galax¬ 
ies except for the youngest stars of age less than 3 x 10 s yr. At first 
sight, it seems surprising that the deviation between the two popula¬ 
tions can have occurred so recently. One interpretation is that galax¬ 
ies in this mass range may transition back and forth between star 
forming and quiescent (i.e. show bursts of star formation) on time- 
scales of this order, and that might also help to explain why the SFR 
of the bluest galaxies seems to increase to the present day. However, 
this seems at odds with the observation that red and blue galaxies 
are observed to have very different metallicities: the stellar-mass- 
weighted metallicity of the young (age less than 2.5 Gyr) stars is 
0.036 in red galaxies, 8 and 0.019 in blue galaxies, which would 
suggest that the two form distinct populations. 

Whatever the interpretation of the observations, it seems un¬ 
likely that they can be made compatible with the model galaxies, 
which show widely divergent SFHs for red and blue galaxies for 
stars younger than 5 Gyr. It would seem that termination of star for- 

8 In units of the mass fraction of metals with respect to Hydrogen; in these 
units solar metallicity is Zq = 0.02 
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Figure 5. The averaged SFH from VESPA (solid lines), the GWB11 mod¬ 
els (dashed-dotted lines), and the HWT14 model (dashed lines) for red and 
blue galaxies with masses greater than 10 11 M.gand within redshift inter¬ 
val 0.10 < 2 < 0.12. 

mation is too abrupt in the models as compared the observations, 
and lacks the possibility of retriggering of star formation at later 
times. 

The comparison of observed and model galaxies is compli¬ 
cated by the effects of metallicity, dust attenuation and finite fi¬ 
bre aperture on the measured colours. In the models, there is a 
very strong correlation between stellar mass and metallicity, and 
between SFR and extinction; in addition; there is no aperture cor¬ 
rection. In real galaxies, the scatter is observed to be much higher, 
and there will be a redshift-dependent colour correction for the fi¬ 
nite aperture. This, together with observational error, is likely to 
move the observed blue and red populations towards each other, so 
caution should be exercised before drawing definitive conclusions. 
On the other hand, SA models struggle to match even the colour 
distribution of galaxies (see, for example. Fig. 9 of HWT14), and 
the distinction between the models and the observed SDSS data is 
so large that it is hard to dismiss it lightly. This issue will be investi¬ 
gated in a subsequent paper, and highlights the power in comparing 
fully-resolved SFHs between models and data. 


5 THE EVOLUTION OF SFHS 

This section looks at the predicted SFHs of galaxies at different 
redshifts. We are interested both in the history of the mean (and me¬ 
dian) galaxy population and of the scatter about that mean. This can 
have important implications for the interpretation of high-redshift 
galaxies that often rely upon postulated SFHs (Boquien et al. 2014; 
Pacifici et al. 2015). 

5.1 Mean SFHs 

Fig. 6 shows the mean, mass-weighted SFRs of all galaxies with 
mass greater than 10 9,s Mq at four different redshifts: 0, 1, 3 and 
7. (Note that this plot differs from previous ones in that we are plot¬ 
ting Af/At rather than Af/A log t; we do this to make it easier to 
detect any decrease in the SFR at recent times). The SFR increases 
rapidly at early times, then slows down, with a decline to late times 
(i.e. low stellar ages) being apparent for 2 < 2. 



Figure 6. The average star formation rates of model galaxies from the 
HWT14 model, with mass greater than 10 9,5 Mq, at four different red¬ 
shifts, as shown. The symbols show the model predictions; the curves are 
fits to the data as described in the text: solid lines, gamma model (Equa¬ 
tion 1); dashed lines, two-power model (Equation 2). 

To characterize the SFR at high redshift, we fit a gamma model 
featuring a power-law increase in star formation at early times fol¬ 
lowed by an exponential decline: 9 

^-=Ax p e~ x , x = (i) 

df r 

Here i a is the age of the galaxy, p sets the rate at which star forma¬ 
tion builds up, and r is the characteristic time-scale over which star 
formation declines. 

At all redshifts above 2 = 1, the SFHs are well-fit by a single 
set of parameters, p = 1.5 and r = 2.0 Gyr, with only f a varying to 
reflect the age of the galaxy. 10 That this is the case is not surprising 
but reflects that fact that the majority of stars in the Universe are 
bom within galaxies whose mass exceeds 3 x 10 9 Mq: each of the 
SFHs shown in Fig. 6 then mirrors the cosmic SFH. Using this set 
of parameters then star formation begins in our model at 2 « 12, 
0.4 Gyr after the big bang, and levels off (i.e. d 2 //df 2 = 0) 3 Gyr 
later, at z ~ 2. 

At lower redshifts, it becomes apparent that an exponential 
decline is too steep. Instead, a two-power model is preferred: 

fL = _d_ x= t -^±. (2) 

dt x~ p + x q ’ r 

Taking p = 1.39, q = 1.33 and r = 3.0 Gyr gives a good fit at all 
redshifts. 

Behroozi et al. (2013) found identical fitting formulae to 
those used here to be good fits to the SFHs of galaxies in their 
abundance matching method to populate haloes with galaxies that 
match observed stellar mass functions and SFRs. However, for 
10 12 Mghaloes, they find a value of q, that determines the rate of 
decay of the SFR at late times, to be significantly higher than that 
quoted above: the reason for this difference is not clear. 

We stress that the curves shown in Fig. 6 are for the mean star 

9 The same functional form with p = 1 was shown by Simha et al. (2014) 
to be a good fit to the individual SFHs of most of their galaxies in SPH 
simulations of galaxy formation. 

10 For 2 > 2 these parameters are degenerate, but we choose to freeze 
them at the values found at lower redshift. 
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Figure 7. The sSFR at different redshifts for the HWT14 model. The sub¬ 
script “first 10 %” refers to an average over the most recent 10 per cent of 
the age of the Universe at that time. The diamonds show the median values 
and the squares show the inverse age of the Universe at each redshift. 

formation rate averaged over all galaxies with mass greater than 
10 9 ' 5 Mq. As is apparent from Fig. 4, high mass galaxies form 
their stars earlier, and low mass galaxies later, than the mean trend. 
We show in Section 5.3 that this is driven primarily by a mass- 
dependence in the cessation of star formation, and in Section 5.4 
that there is considerable variation between individual galaxies. 

5.2 Specific SFRs and quiescent fractions 

Fig. 7 shows the specific SFR (sSFR) of all galaxies with M > 
3 x 10 9 M 0 in the F1WT14 mode at several different redshifts. 
The SFR here is averaged over the most recent 10 per cent of the 
age of the Universe at that time. The diamond symbols show the 
location of the median values, and the squares show the inverse age 
of the Universe at that redshift. 

There is a small spread around the modal sSFR, with the vast 
majority of star forming galaxies lying within about ±0.25 dex of 
the peak. Flowever, there is a long tail of galaxies extending to low 
sSFRs, which becomes more prominent at low redshifts; indeed, at 
z = 0, 35 per cent lie off the left-hand edge of the plot altogether. 
That is why, below z = 1, the median values lie well to the left of 
the mode. 

Above 2 « 0.35, most galaxies are forming stars at a rate that 
would more than double their mass in the age of the Universe; 11 
below that redshift, the opposite is true. It is at this time that there 
is a strong shift from star forming to non-star forming galaxies. 
There is no sharp distinction between the two but, using Fig. 7 as 
a guide, we define a galaxy to be quiescent if it has formed fewer 
than 3 per cent of its stars in the most recent 10 per cent of the age 
of the Universe, t z , at that redshift, i.e. sSFRg rst 10 % < 0.3/t z . 

Fig. 8 shows the fraction of quiescent galaxies as a function 
of mass and redshift. At high redshift, there is some suppression 
of star formation in dwarf galaxies, reflecting the strong feedback 
from supernova in the model. However, only a small fraction of 
galaxies are affected and there is a much greater growth of the pas¬ 
sive population at redshifts below z ~ 3. Once again, there is clear 

11 Note that this is a number-weighted average, so that does not mean that 
mean SFR peaked at that time. 



Figure 8 . The quiescent fraction of galaxies in the F1WT14 model, in several 
different mass bins, and at several different redshifts, as shown. 


Table 3. The median sSFRs, and the quiescent fraction of galaxies in the 
HWT14 model with mass exceeding 3 X 10 9 Mq. The columns are: red¬ 
shift; age of the Universe in Gyr, t/Gyr; median specific SFR multiplied 
by the age of the Universe, So.5 all! th e same but restricted to star form¬ 
ing galaxies, So .5 s f; quiescent fraction weighted by number, Q num ; and 
quiescent fraction weighted by mass, Q m ass ■ 


Redshift 

t/Gyr 

■So. 5 , all 

*S'o.5,sf 

Qnum 

Qmass 

6.97 

0.77 

3.71 

3.71 

0.00 

0.00 

5.03 

1.16 

3.17 

3.17 

0.00 

0.00 

3.95 

1.56 

2.81 

2.82 

0.01 

0.01 

3.11 

2.06 

2.50 

2.53 

0.03 

0.02 

2.07 

3.17 

2.02 

2.09 

0.07 

0.08 

1.04 

5.69 

1.43 

1.60 

0.18 

0.28 

0.35 

9.79 

0.99 

1.31 

0.33 

0.48 

0.00 

13.80 

0.69 

1.17 

0.41 

0.58 


evidence of down-sizing in that more massive galaxies start to be¬ 
come quiescent earlier than lower-mass ones. 

Table 3 lists several measures of star formation activity. 
So. 5 ,ail is the median sSFR of the sample multiplied by the age 
of the Universe at that redshift; likewise So.s.af is the same thing, 
but restricted to star forming galaxies. The two begin to differ sig¬ 
nificantly below a redshift of about 3 once the quiescent fraction 
begins to rise. This fraction is listed weighted both by galaxy num¬ 
ber and by galaxy mass, from which it can be seen that about 58 
per cent of stars in the current-day Universe lie in galaxies that are 
not actively star forming. 

The medians listed in the table show that there has been a 
steady decline in star formation activity in galaxies from a redshift 
of at least 7 right through to the current day. However, even as re¬ 
cently as 2 « 0.35, most galaxies were still forming stars at a rate 
that would more than double their mass within the age of the Uni¬ 
verse at that time. 


5.3 The cause of mass-dependent SFHs 

The top panel of Fig. 9 shows the SFHs of galaxies at 2 = 0 split 
into four different mass bins. In order to better illustrate the onset 
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Figure 9. The SFR as a function of time for the HWT14 model at 2 = 0, 
split into 4 different mass bins, as shown. The upper panel shows all galaxies 
in each mass-range; the middle panel shows star forming galaxies; and the 
lower panel panel shows quiescent galaxies, as defined in Section 5.2. 

of star formation, which is poorly resolved at 2 = 0 using our 
default number of bins, for this section only we use iV max = 4 
(see Section 2.3), giving 29 time-bins at 2 = 0. Age downsizing 
is clearly visible with more massive galaxies forming their stars 
earlier than lower-mass ones. 

To better understand the cause of this, the middle and lower 
panels show the same curves for star forming, and for quiescent 
galaxies, respectively. Although the correspondence is not perfect, 
the agreement between the different mass-bins is much tighter in 
the central panel than in the upper one. That strongly suggests that 
the SFHs of star forming galaxies are very similar, independent 
of the mass of the galaxy, and that the primary driver of the mass- 
dependence is the different evolution of the quiescent fraction. Note 
that we tried only a single definition of quiescence and it is likely 
that the residual mass-dependence in the central panel could be re¬ 
duced even further if we optimised the definition for that purpose. 

5.4 Individual SFHs 

In this section, we restrict our attention to star forming galaxies. 

Although the mean SFH is well-described by a simple func¬ 
tional form, Fig. 10 shows that individual galaxies have a wide vari¬ 
ety of histories. This figure shows histograms of the rate of decline 
of the SFR measured by the ratio in two successive time-bins: the 
most recent 10 per cent, and the next most recent 10 per cent, of the 
age of the Universe at that redshift. Galaxies to the left/right of the 
vertical line have declining/increasing SFRs, respectively. 

The distribution of ratios is shown in Fig. 10 for a variety of 



'°9lo( S ^^^tirst 10x/ S S^second 10*) 

Figure 10. The ratio of the SFR in the most recent (’first’) 10 per cent, to 
that in the previous (’second’) 10 per cent, of the age of the Universe at 
that redshift, for galaxies with masses exceeding 10 9 ' 5 Mq in the HWT14 
model. The j/-scale corresponds to the 2 = 0 curve: the higher-redshift 
curves are offset by successive factors of two to space them out in the in¬ 
direction. The diamond symbols show the median values. 


redshifts. There is a gradual shift from increasing to decreasing spe¬ 
cific star formation rates as the Universe ages. When measured in 
this way, an equal balance between increasing and decreasing SFRs 
is achieved somewhere between redshifts 1 and 2. At all times, 
however, there is a significant fraction of galaxies lying in each 
of these populations: at 2 = 7, four-fifths of galaxies show an in¬ 
creasing SFR. and at 2 = 0 three-quarters show a decreasing one. 


6 CONCLUSIONS 

In this paper, we have introduced the recording of star formation 
histories (SFHs) in the L-GALAXIES SA model. At any given 
point in a galaxy’s evolution, the mass of recently-formed stars is 
recorded in bins of time resolution equal to that of the timestep in 
the SA model (l-2xl0 7 yr). These bins are gradually merged to¬ 
gether as the galaxy ages, such that older stars are grouped together 
into larger bins. 

We investigate the extent to which SFHs may be used to re¬ 
construct stellar spectra in post-processing; we compare our SFHs 
to those in the publicly-available VESPA catalogue extracted from 
SDSS-DR7 data; and we investigate in our favoured SA model 
(Henriques et al. 2014) the evolution of SFHs as a function of 
galaxy mass. Our key results are as follows: 

• Post-processing reconstruction of magnitudes in various ob¬ 
servational bands gives good agreement with on-the-fly accumula¬ 
tion of luminosity, provided that JVmax sS 2 (which equates to 16 
bins at 2 = 0). Quantitatively, the rms difference between raw and 
reconstructed magnitudes is less than 0.05 for all bands except the 
far-UV, for which it is 0.29. 

• The SA models show reasonable qualitative agreement with 
the observed SFHs of the SDSS Main Galaxy Sample from the 
VESPA catalogue, with the GWB11 model fitting better for stars 
older than 1 Gyr, and the HWT14 model fitting better for younger 
stars. 

• When divided the mass, both the observations and models 


© 2013 RAS, MNRAS 000, 1-11 


















10 Shamshiri et al. 


show a trend for more massive galaxies to form their stars earlier 
and have lower current sSFRs than lower-mass galaxies. 

• When divided by colour, the agreement is poorer. Both ver¬ 
sions of the SA model show much more extreme variation in SFH 
with colour than do observed galaxies from the VESPA data base. 
In the model, the SFHs of red and blue galaxies begin to differ as 
long ago as 5 Gyr, compared to just 0.3 Gyr for observed galaxies. 
One possible explanation could be that real galaxies show repeated 
episodes of star formation that are not present in the models. We 
note, however, that a more rigorous investigation of the data is re¬ 
quired before drawing any definitive conclusions. 

• At 2 ^ 1 the mean SFR of all model galaxies with stellar 
mass greater than 3 x 10 9 Mgis well fitted by the formula df/dt oc 
x 1 '°e ~ x , where x = (f a — £)/2.0 Gyr. Here t is the lookback time 
and f a is the age of the galaxy. At later times, the SFR declines less 
rapidly and a two-power model (that contains an extra parameter) is 
a better fit over the whole of cosmic history: d//d£ oc 1/(a; -1 ' 39 + 
a; 1 ’ 33 ), where a: = (£ a — £)/3.0Gyr. 

• Although star formation rates have been declining for more 
than half the history of the Universe, the typical (median) star form¬ 
ing galaxy today is still forming stars at a rate that will more than 
double its mass in a Hubble time. 

• We define a galaxy to be quiescent if it forms fewer than 3 
per cent of its stars in the most recent 10 per cent of the age of the 
Universe, t z at that redshift, i.e. sSFR< 0.3 /t z . Then the quiescent 
fraction begins to increase rapidly below z ~ 3, reaching 41 % by 
number and 58 % by mass for galaxies with M > 10 9 Mq at the 
current day. 

• Our model produces a small fraction of quiescent dwarf galax¬ 
ies, M < 10 10 Mgat all times, but the main effect that we see is 
consistent with downsizing in that more massive galaxies become 
quiescent first, followed by successively lower mass galaxies as the 
Universe ages. 

• When split by mass, and using a finer time resolution in the 
SFHs, downsizing is very clear to see, with the peak of the SFR 
shifting from a lookback time of about 11 Gyr in the most mas¬ 
sive galaxies (current-day mass greater than 10 11 Mq) to less than 
8 Gyr in lower-mass systems (3 x 10 9 -10 10 Mq). 

• When split into quiescent and a star forming populations, the 
differences between the mean SFHs of star forming galaxies of dif¬ 
ferent mass is much reduced. Downsizing thus has its origin in an 
earlier transition from star forming to quiescent status in galaxies 
that are more massive at the current-day. 

• Although the mean SFHs are well defined, there is a huge dis¬ 
persion in the SFHs of individual galaxies such that, even at the 
current day, many galaxies still have increasing SFRs. At z = 7, 
four-fifths of galaxies show an increasing SFR, and at z = 0 three- 
quarters show a decreasing one; an equal balance between galaxies 
with increasing and decreasing SFRs is achieved somewhere be¬ 
tween redshifts 1 and 2. 

As can be seen from the above, one of the key drivers of galaxy 
evolution is the rate at which star formation is quenched. This is 
investigated in a companion paper, Henriques et al. (in prepara¬ 
tion), that undertakes a detailed comparison with observations of 
the quenched fraction as a function of environment and mass. There 
it is shown that the HWT14 model does a much better job than pre¬ 
vious incarnations of the L-GALAXIES SA model in terminating 
star formation in massive galaxies, whilst allowing continued star 
formation in low-mass satellites, though the quantitative agreement 
is still far from ideal. 

An earlier paper, Yates et al. (2013), combined the SFHs with 


a multicomponent model for stellar feedback to investigate the 
metallicity evolution of galaxies. This then enables us to construct 
metallicity histories for galaxies along the lines of the SFHs pre¬ 
sented in this paper. Unfortunately, the observational data from 
VESPA is currently unable to constrain the metallicity histories 
with any degree of certainty. 

The low-resolution (JV max = 2) SFHs for the HWT14 SA 
model, presented in this paper, are publicly available to download 
from the Millennium data base 12 and have been used to reconstruct 
predicted fluxes in post-processing. Higher resolution catalogues 
are available from the authors upon request. 
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APPENDIX A: THE CHOICE OF VESPA CATALOGUE 

The VESPA SDSS-7 catalogue (Tojeiro et al. 2009) contains galax¬ 
ies with a wide variety of data quality, some showing reconstructed 
mass errors that are greater than 100 %. In order not to bias the 
results, we include the whole sample in our analysis. We show av¬ 
erage SFHs weighted by galaxy number, rather than galaxy mass, 
so as to minimize the effect of the errors in the mass reconstruc¬ 
tion. We have checked that restricting the analysis to the galaxies 
with the best data quality does, in fact, lead to qualitatively similar 
SFHs. 

The average SFHs of galaxies in this subsample, weighted 
by galaxy number, are shown in Fig. Al for two different SEDs 
(Bruzual & Chariot 2003, hereafter BC03 and Maraston 2005, here- 
afterM05) and two different dust models. The one-component dust 
model is a uniform screen applied to the whole stellar popula¬ 
tion; the two-component model adds in extra absorption in front 
of young stars. 

First note that the M05 models show a much smoother change 
in the SFR between look-back times of 0. l-10Gyr than do those 
of BC03. In such a large galaxy sample, it is hard to think of a 



Figure Al. The average SFH of VESPA galaxies obtained using the SEDs 
of BC03 (solid lines) and M05 (dashed lines). Blue lines show results for a 
one-parameter dust model whereas red shows a two-parameter dust model. 
The solid, black line has a slope of unity, corresponding to a constant SFR. 

plausible reason for this and the SA models show no such feature. 
For that reason, we use the M05 results. 

Both the BC03 and, to a lesser extent, the M05 results for the 
two-dust model show a significant increase in SFR at ages less than 
0.1 Gyr. Again, this seems implausible and suggests that there is not 
enough constraining power in the data: the model has presumably 
confused dust obscured, young stars with some older population, 
perhaps to explain some spectral feature that is not well fitted by 
the SEDs. 

Throughout the body of the paper, we use the VESPA results 
for M05 and a single dust model. 
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